/*
 * gradient_descent.h
 *
 *  Created on: Jun 3, 2011
 *      Author: tqlong
 */

#ifndef GRADIENT_DESCENT_H_
#define GRADIENT_DESCENT_H_

#include "line_search.h"

template<typename Function>
  class GradientDescent
  {
  public:
    typedef Function func_type;
    typedef LineSearch<Function> linesearch_type;
    typedef typename Function::range_type range_type;
    typedef typename Function::domain_type domain_type;
  public:
    GradientDescent(func_type &f, const domain_type &x0, domain_type &x) :
      f_(f), x0_(x0), x_(&x), step_(1.0)
    {
      if (!f.has_value_)
        {
          std::cerr << "Function needs to produce value() at " << __FILE__
              << ":" << __LINE__ << std::endl;
          exit(1);
        }
      if (!f.has_gradient_)
        {
          std::cerr << "Function needs to produce gradient() at " << __FILE__
              << ":" << __LINE__ << std::endl;
          exit(1);
        }
    }

    bool
    search()
    {
      *x_ = x0_;
      fx_ = f_.value(*x_);
      iter_ = 0;
      for (iter_ = 0; iter_ < 100; iter_++)
        {
          grad_ = f_.gradient(*x_);
          domain_type newx, dir = -grad_;
          linesearch_type ls(f_, *x_, newx, fx_, grad_, dir, step_);

          bool ok = ls.search();
          *x_ = newx;
          fx_ = ls.fx();
          step_ = ls.step();

          double grad_norm = arma::norm(grad_, 2);

          std::cout << "iter=" << iter_ << " fx=" << fx_ << " step=" << step_
              << " |grad|=" << grad_norm << std::endl;
          if (!ok || grad_norm < 1e-5) break; // line search failed
        }
      return true;
    }

    void
    value_gradient(const domain_type& x, range_type& value, domain_type& grad)
    {
      if (f_.has_both_)
        f_.value_gradient(x, value, grad);
      else
        {
          value = f_.value(x);
          grad = f_.gradient(x);
        }
    }

    range_type
    fx() const
    {
      return fx_;
    }

    const domain_type&
    grad() const
    {
      return grad_;
    }
  protected:
    func_type &f_;
    const domain_type &x0_;
    domain_type *x_;
    range_type fx_;
    domain_type grad_;
    double step_;
    int iter_;
  };

#endif /* GRADIENT_DESCENT_H_ */
